home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / LINUX / MATH_EMU.ZIP / MATH_EMU / FPU_LOG.C < prev    next >
Encoding:
C/C++ Source or Header  |  1979-12-31  |  21.4 KB  |  599 lines
  1. /*        $NetBSD$  */
  2.  
  3. /*
  4.  * Copyright (c) 1995  Ken Nakata
  5.  *        All rights reserved.
  6.  *
  7.  * Redistribution and use in source and binary forms, with or without
  8.  * modification, are permitted provided that the following conditions
  9.  * are met:
  10.  * 1. Redistributions of source code must retain the above copyright
  11.  *    notice, this list of conditions and the following disclaimer.
  12.  * 2. Redistributions in binary form must reproduce the above copyright
  13.  *    notice, this list of conditions and the following disclaimer in the
  14.  *    documentation and/or other materials provided with the distribution.
  15.  * 3. Neither the name of the author nor the names of its contributors
  16.  *    may be used to endorse or promote products derived from this software
  17.  *    without specific prior written permission.
  18.  *
  19.  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  20.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  21.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  22.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  23.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  24.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  25.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  26.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  27.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  28.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  29.  * SUCH DAMAGE.
  30.  *
  31.  *        @(#)fpu_log.c       10/8/95
  32.  */
  33.  
  34. #include "types.h"
  35.  
  36. #include "fpu_emul.h"
  37.  
  38. static u_int logA6[] = { 0x3FC2499A, 0xB5E4040B };
  39. static u_int logA5[] = { 0xBFC555B5, 0x848CB7DB };
  40. static u_int logA4[] = { 0x3FC99999, 0x987D8730 };
  41. static u_int logA3[] = { 0xBFCFFFFF, 0xFF6F7E97 };
  42. static u_int logA2[] = { 0x3FD55555, 0x555555A4 };
  43. static u_int logA1[] = { 0xBFE00000, 0x00000008 };
  44.  
  45. static u_int logB5[] = { 0x3F175496, 0xADD7DAD6 };
  46. static u_int logB4[] = { 0x3F3C71C2, 0xFE80C7E0 };
  47. static u_int logB3[] = { 0x3F624924, 0x928BCCFF };
  48. static u_int logB2[] = { 0x3F899999, 0x999995EC };
  49. static u_int logB1[] = { 0x3FB55555, 0x55555555 };
  50.  
  51. /* sfpn = shortened fp number; can represent only positive numbers */
  52. static struct sfpn {
  53.     int             sp_exp;
  54.     u_int sp_m0, sp_m1;
  55. } logtbl[] = {
  56.     { 0x3FFE - 0x3fff, 0xFE03F80FU, 0xE03F80FEU },
  57.     { 0x3FF7 - 0x3fff, 0xFF015358U, 0x833C47E2U },
  58.     { 0x3FFE - 0x3fff, 0xFA232CF2U, 0x52138AC0U },
  59.     { 0x3FF9 - 0x3fff, 0xBDC8D83EU, 0xAD88D549U },
  60.     { 0x3FFE - 0x3fff, 0xF6603D98U, 0x0F6603DAU },
  61.     { 0x3FFA - 0x3fff, 0x9CF43DCFU, 0xF5EAFD48U },
  62.     { 0x3FFE - 0x3fff, 0xF2B9D648U, 0x0F2B9D65U },
  63.     { 0x3FFA - 0x3fff, 0xDA16EB88U, 0xCB8DF614U },
  64.     { 0x3FFE - 0x3fff, 0xEF2EB71FU, 0xC4345238U },
  65.     { 0x3FFB - 0x3fff, 0x8B29B775U, 0x1BD70743U },
  66.     { 0x3FFE - 0x3fff, 0xEBBDB2A5U, 0xC1619C8CU },
  67.     { 0x3FFB - 0x3fff, 0xA8D839F8U, 0x30C1FB49U },
  68.     { 0x3FFE - 0x3fff, 0xE865AC7BU, 0x7603A197U },
  69.     { 0x3FFB - 0x3fff, 0xC61A2EB1U, 0x8CD907ADU },
  70.     { 0x3FFE - 0x3fff, 0xE525982AU, 0xF70C880EU },
  71.     { 0x3FFB - 0x3fff, 0xE2F2A47AU, 0xDE3A18AFU },
  72.     { 0x3FFE - 0x3fff, 0xE1FC780EU, 0x1FC780E2U },
  73.     { 0x3FFB - 0x3fff, 0xFF64898EU, 0xDF55D551U },
  74.     { 0x3FFE - 0x3fff, 0xDEE95C4CU, 0xA037BA57U },
  75.     { 0x3FFC - 0x3fff, 0x8DB956A9U, 0x7B3D0148U },
  76.     { 0x3FFE - 0x3fff, 0xDBEB61EEU, 0xD19C5958U },
  77.     { 0x3FFC - 0x3fff, 0x9B8FE100U, 0xF47BA1DEU },
  78.     { 0x3FFE - 0x3fff, 0xD901B203U, 0x6406C80EU },
  79.     { 0x3FFC - 0x3fff, 0xA9372F1DU, 0x0DA1BD17U },
  80.     { 0x3FFE - 0x3fff, 0xD62B80D6U, 0x2B80D62CU },
  81.     { 0x3FFC - 0x3fff, 0xB6B07F38U, 0xCE90E46BU },
  82.     { 0x3FFE - 0x3fff, 0xD3680D36U, 0x80D3680DU },
  83.     { 0x3FFC - 0x3fff, 0xC3FD0329U, 0x06488481U },
  84.     { 0x3FFE - 0x3fff, 0xD0B69FCBU, 0xD2580D0BU },
  85.     { 0x3FFC - 0x3fff, 0xD11DE0FFU, 0x15AB18CAU },
  86.     { 0x3FFE - 0x3fff, 0xCE168A77U, 0x25080CE1U },
  87.     { 0x3FFC - 0x3fff, 0xDE1433A1U, 0x6C66B150U },
  88.     { 0x3FFE - 0x3fff, 0xCB8727C0U, 0x65C393E0U },
  89.     { 0x3FFC - 0x3fff, 0xEAE10B5AU, 0x7DDC8ADDU },
  90.     { 0x3FFE - 0x3fff, 0xC907DA4EU, 0x871146ADU },
  91.     { 0x3FFC - 0x3fff, 0xF7856E5EU, 0xE2C9B291U },
  92.     { 0x3FFE - 0x3fff, 0xC6980C69U, 0x80C6980CU },
  93.     { 0x3FFD - 0x3fff, 0x82012CA5U, 0xA68206D7U },
  94.     { 0x3FFE - 0x3fff, 0xC4372F85U, 0x5D824CA6U },
  95.     { 0x3FFD - 0x3fff, 0x882C5FCDU, 0x7256A8C5U },
  96.     { 0x3FFE - 0x3fff, 0xC1E4BBD5U, 0x95F6E947U },
  97.     { 0x3FFD - 0x3fff, 0x8E44C60BU, 0x4CCFD7DEU },
  98.     { 0x3FFE - 0x3fff, 0xBFA02FE8U, 0x0BFA02FFU },
  99.     { 0x3FFD - 0x3fff, 0x944AD09EU, 0xF4351AF6U },
  100.     { 0x3FFE - 0x3fff, 0xBD691047U, 0x07661AA3U },
  101.     { 0x3FFD - 0x3fff, 0x9A3EECD4U, 0xC3EAA6B2U },
  102.     { 0x3FFE - 0x3fff, 0xBB3EE721U, 0xA54D880CU },
  103.     { 0x3FFD - 0x3fff, 0xA0218434U, 0x353F1DE8U },
  104.     { 0x3FFE - 0x3fff, 0xB92143FAU, 0x36F5E02EU },
  105.     { 0x3FFD - 0x3fff, 0xA5F2FCABU, 0xBBC506DAU },
  106.     { 0x3FFE - 0x3fff, 0xB70FBB5AU, 0x19BE3659U },
  107.     { 0x3FFD - 0x3fff, 0xABB3B8BAU, 0x2AD362A5U },
  108.     { 0x3FFE - 0x3fff, 0xB509E68AU, 0x9B94821FU },
  109.     { 0x3FFD - 0x3fff, 0xB1641795U, 0xCE3CA97BU },
  110.     { 0x3FFE - 0x3fff, 0xB30F6352U, 0x8917C80BU },
  111.     { 0x3FFD - 0x3fff, 0xB7047551U, 0x5D0F1C61U },
  112.     { 0x3FFE - 0x3fff, 0xB11FD3B8U, 0x0B11FD3CU },
  113.     { 0x3FFD - 0x3fff, 0xBC952AFEU, 0xEA3D13E1U },
  114.     { 0x3FFE - 0x3fff, 0xAF3ADDC6U, 0x80AF3ADEU },
  115.     { 0x3FFD - 0x3fff, 0xC2168ED0U, 0xF458BA4AU },
  116.     { 0x3FFE - 0x3fff, 0xAD602B58U, 0x0AD602B6U },
  117.     { 0x3FFD - 0x3fff, 0xC788F439U, 0xB3163BF1U },
  118.     { 0x3FFE - 0x3fff, 0xAB8F69E2U, 0x8359CD11U },
  119.     { 0x3FFD - 0x3fff, 0xCCECAC08U, 0xBF04565DU },
  120.     { 0x3FFE - 0x3fff, 0xA9C84A47U, 0xA07F5638U },
  121.     { 0x3FFD - 0x3fff, 0xD2420487U, 0x2DD85160U },
  122.     { 0x3FFE - 0x3fff, 0xA80A80A8U, 0x0A80A80BU },
  123.     { 0x3FFD - 0x3fff, 0xD7894992U, 0x3BC3588AU },
  124.     { 0x3FFE - 0x3fff, 0xA655C439U, 0x2D7B73A8U },
  125.     { 0x3FFD - 0x3fff, 0xDCC2C4B4U, 0x9887DACCU },
  126.     { 0x3FFE - 0x3fff, 0xA4A9CF1DU, 0x96833751U },
  127.     { 0x3FFD - 0x3fff, 0xE1EEBD3EU, 0x6D6A6B9EU },
  128.     { 0x3FFE - 0x3fff, 0xA3065E3FU, 0xAE7CD0E0U },
  129.     { 0x3FFD - 0x3fff, 0xE70D785CU, 0x2F9F5BDCU },
  130.     { 0x3FFE - 0x3fff, 0xA16B312EU, 0xA8FC377DU },
  131.     { 0x3FFD - 0x3fff, 0xEC1F392CU, 0x5179F283U },
  132.     { 0x3FFE - 0x3fff, 0x9FD809FDU, 0x809FD80AU },
  133.     { 0x3FFD - 0x3fff, 0xF12440D3U, 0xE36130E6U },
  134.     { 0x3FFE - 0x3fff, 0x9E4CAD23U, 0xDD5F3A20U },
  135.     { 0x3FFD - 0x3fff, 0xF61CCE92U, 0x346600BBU },
  136.     { 0x3FFE - 0x3fff, 0x9CC8E160U, 0xC3FB19B9U },
  137.     { 0x3FFD - 0x3fff, 0xFB091FD3U, 0x8145630AU },
  138.     { 0x3FFE - 0x3fff, 0x9B4C6F9EU, 0xF03A3CAAU },
  139.     { 0x3FFD - 0x3fff, 0xFFE97042U, 0xBFA4C2ADU },
  140.     { 0x3FFE - 0x3fff, 0x99D722DAU, 0xBDE58F06U },
  141.     { 0x3FFE - 0x3fff, 0x825EFCEDU, 0x49369330U },
  142.     { 0x3FFE - 0x3fff, 0x9868C809U, 0x868C8098U },
  143.     { 0x3FFE - 0x3fff, 0x84C37A7AU, 0xB9A905C9U },
  144.     { 0x3FFE - 0x3fff, 0x97012E02U, 0x5C04B809U },
  145.     { 0x3FFE - 0x3fff, 0x87224C2EU, 0x8E645FB7U },
  146.     { 0x3FFE - 0x3fff, 0x95A02568U, 0x095A0257U },
  147.     { 0x3FFE - 0x3fff, 0x897B8CACU, 0x9F7DE298U },
  148.     { 0x3FFE - 0x3fff, 0x94458094U, 0x45809446U },
  149.     { 0x3FFE - 0x3fff, 0x8BCF55DEU, 0xC4CD05FEU },
  150.     { 0x3FFE - 0x3fff, 0x92F11384U, 0x0497889CU },
  151.     { 0x3FFE - 0x3fff, 0x8E1DC0FBU, 0x89E125E5U },
  152.     { 0x3FFE - 0x3fff, 0x91A2B3C4U, 0xD5E6F809U },
  153.     { 0x3FFE - 0x3fff, 0x9066E68CU, 0x955B6C9BU },
  154.     { 0x3FFE - 0x3fff, 0x905A3863U, 0x3E06C43BU },
  155.     { 0x3FFE - 0x3fff, 0x92AADE74U, 0xC7BE59E0U },
  156.     { 0x3FFE - 0x3fff, 0x8F1779D9U, 0xFDC3A219U },
  157.     { 0x3FFE - 0x3fff, 0x94E9BFF6U, 0x15845643U },
  158.     { 0x3FFE - 0x3fff, 0x8DDA5202U, 0x37694809U },
  159.     { 0x3FFE - 0x3fff, 0x9723A1B7U, 0x20134203U },
  160.     { 0x3FFE - 0x3fff, 0x8CA29C04U, 0x6514E023U },
  161.     { 0x3FFE - 0x3fff, 0x995899C8U, 0x90EB8990U },
  162.     { 0x3FFE - 0x3fff, 0x8B70344AU, 0x139BC75AU },
  163.     { 0x3FFE - 0x3fff, 0x9B88BDAAU, 0x3A3DAE2FU },
  164.     { 0x3FFE - 0x3fff, 0x8A42F870U, 0x5669DB46U },
  165.     { 0x3FFE - 0x3fff, 0x9DB4224FU, 0xFFE1157CU },
  166.     { 0x3FFE - 0x3fff, 0x891AC73AU, 0xE9819B50U },
  167.     { 0x3FFE - 0x3fff, 0x9FDADC26U, 0x8B7A12DAU },
  168.     { 0x3FFE - 0x3fff, 0x87F78087U, 0xF78087F8U },
  169.     { 0x3FFE - 0x3fff, 0xA1FCFF17U, 0xCE733BD4U },
  170.     { 0x3FFE - 0x3fff, 0x86D90544U, 0x7A34ACC6U },
  171.     { 0x3FFE - 0x3fff, 0xA41A9E8FU, 0x5446FB9FU },
  172.     { 0x3FFE - 0x3fff, 0x85BF3761U, 0x2CEE3C9BU },
  173.     { 0x3FFE - 0x3fff, 0xA633CD7EU, 0x6771CD8BU },
  174.     { 0x3FFE - 0x3fff, 0x84A9F9C8U, 0x084A9F9DU },
  175.     { 0x3FFE - 0x3fff, 0xA8489E60U, 0x0B435A5EU },
  176.     { 0x3FFE - 0x3fff, 0x83993052U, 0x3FBE3368U },
  177.     { 0x3FFE - 0x3fff, 0xAA59233CU, 0xCCA4BD49U },
  178.     { 0x3FFE - 0x3fff, 0x828CBFBEU, 0xB9A020A3U },
  179.     { 0x3FFE - 0x3fff, 0xAC656DAEU, 0x6BCC4985U },
  180.     { 0x3FFE - 0x3fff, 0x81848DA8U, 0xFAF0D277U },
  181.     { 0x3FFE - 0x3fff, 0xAE6D8EE3U, 0x60BB2468U },
  182.     { 0x3FFE - 0x3fff, 0x80808080U, 0x80808081U },
  183.     { 0x3FFE - 0x3fff, 0xB07197A2U, 0x3C46C654U },
  184. };
  185.  
  186. static struct fpn *__fpu_logn __P((struct fpemu *fe));
  187.  
  188. /*
  189.  * natural log - algorithm taken from Motorola FPSP,
  190.  * except this doesn't bother to check for invalid input.
  191.  */
  192. static struct fpn *
  193. __fpu_logn(fe)
  194.      struct fpemu *fe;
  195. {
  196.     static struct fpn X, F, U, V, W, KLOG2;
  197.     struct fpn *d;
  198.     int i, k;
  199.  
  200.     CPYFPN(&X, &fe->fe_f2);
  201.  
  202.     /* see if |X-1| < 1/16 approx. */
  203.     if ((-1 == X.fp_exp && (0xf07d0000U >> (31 - FP_LG)) <= X.fp_mant[0]) ||
  204.           (0 == X.fp_exp && X.fp_mant[0] <= (0x88410000U >> (31 - FP_LG)))) {
  205.           /* log near 1 */
  206. #ifdef DEBUG
  207.           if (fpu_debug_level & DL_ARITH)
  208.               printf("__fpu_logn: log near 1\n");
  209. #endif
  210.           fpu_const(&fe->fe_f1, 0x32);
  211.           /* X+1 */
  212.           d = fpu_add(fe);
  213.           CPYFPN(&V, d);
  214.  
  215.           CPYFPN(&fe->fe_f1, &X);
  216.           fpu_const(&fe->fe_f2, 0x32); /* 1.0 */
  217.           fe->fe_f2.fp_sign = 1; /* -1.0 */
  218.           /* X-1 */
  219.           d = fpu_add(fe);
  220.           CPYFPN(&fe->fe_f1, d);
  221.           /* 2(X-1) */
  222.           fe->fe_f1.fp_exp++; /* *= 2 */
  223.           CPYFPN(&fe->fe_f2, &V);
  224.           /* U=2(X-1)/(X+1) */
  225.           d = fpu_div(fe);
  226.           CPYFPN(&U, d);
  227.           CPYFPN(&fe->fe_f1, d);
  228.           CPYFPN(&fe->fe_f2, d);
  229.           /* V=U*U */
  230.           d = fpu_mul(fe);
  231.           CPYFPN(&V, d);
  232.           CPYFPN(&fe->fe_f1, d);
  233.           CPYFPN(&fe->fe_f2, d);
  234.           /* W=V*V */
  235.           d = fpu_mul(fe);
  236.           CPYFPN(&W, d);
  237.  
  238.           /* calculate U+U*V*([B1+W*(B3+W*B5)]+[V*(B2+W*B4)]) */
  239.  
  240.           /* B1+W*(B3+W*B5) part */
  241.           CPYFPN(&fe->fe_f1, d);
  242.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB5);
  243.           /* W*B5 */
  244.           d = fpu_mul(fe);
  245.           CPYFPN(&fe->fe_f1, d);
  246.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB3);
  247.           /* B3+W*B5 */
  248.           d = fpu_add(fe);
  249.           CPYFPN(&fe->fe_f1, d);
  250.           CPYFPN(&fe->fe_f2, &W);
  251.           /* W*(B3+W*B5) */
  252.           d = fpu_mul(fe);
  253.           CPYFPN(&fe->fe_f1, d);
  254.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB1);
  255.           /* B1+W*(B3+W*B5) */
  256.           d = fpu_add(fe);
  257.           CPYFPN(&X, d);
  258.  
  259.           /* [V*(B2+W*B4)] part */
  260.           CPYFPN(&fe->fe_f1, &W);
  261.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB4);
  262.           /* W*B4 */
  263.           d = fpu_mul(fe);
  264.           CPYFPN(&fe->fe_f1, d);
  265.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB2);
  266.           /* B2+W*B4 */
  267.           d = fpu_add(fe);
  268.           CPYFPN(&fe->fe_f1, d);
  269.           CPYFPN(&fe->fe_f2, &V);
  270.           /* V*(B2+W*B4) */
  271.           d = fpu_mul(fe);
  272.           CPYFPN(&fe->fe_f1, d);
  273.           CPYFPN(&fe->fe_f2, &X);
  274.           /* B1+W*(B3+W*B5)+V*(B2+W*B4) */
  275.           d = fpu_add(fe);
  276.           CPYFPN(&fe->fe_f1, d);
  277.           CPYFPN(&fe->fe_f2, &V);
  278.           /* V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
  279.           d = fpu_mul(fe);
  280.           CPYFPN(&fe->fe_f1, d);
  281.           CPYFPN(&fe->fe_f2, &U);
  282.           /* U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
  283.           d = fpu_mul(fe);
  284.           CPYFPN(&fe->fe_f1, d);
  285.           CPYFPN(&fe->fe_f2, &U);
  286.           /* U+U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
  287.           d = fpu_add(fe);
  288.     } else /* the usual case */ {
  289. #ifdef DEBUG
  290.           if (fpu_debug_level & DL_ARITH)
  291.               printf("__fpu_logn: the usual case. X=(%d,%08x,%08x...)\n",
  292.                        X.fp_exp, X.fp_mant[0], X.fp_mant[1]);
  293. #endif
  294.           k = X.fp_exp;
  295.           /* X <- Y */
  296.           X.fp_exp = fe->fe_f2.fp_exp = 0;
  297.  
  298.           /* get the most significant 7 bits of X */
  299.           F.fp_class = FPC_NUM;
  300.           F.fp_sign = 0;
  301.           F.fp_exp = X.fp_exp;
  302.           F.fp_mant[0] = X.fp_mant[0] & (0xfe000000U >> (31 - FP_LG));
  303.           F.fp_mant[0] |= (0x01000000U >> (31 - FP_LG));
  304.           F.fp_mant[1] = F.fp_mant[2] = F.fp_mant[3] = 0;
  305.           F.fp_sticky = 0;
  306. #ifdef DEBUG
  307.           if (fpu_debug_level & DL_ARITH) {
  308.               printf("__fpu_logn: X=Y*2^k=(%d,%08x,%08x...)*2^%d\n",
  309.                        fe->fe_f2.fp_exp, fe->fe_f2.fp_mant[0],
  310.                        fe->fe_f2.fp_mant[1], k);
  311.               printf("__fpu_logn: F=(%d,%08x,%08x...)\n",
  312.                        F.fp_exp, F.fp_mant[0], F.fp_mant[1]);
  313.           }
  314. #endif
  315.           /* index to the table */
  316.           i = (F.fp_mant[0] >> (FP_LG - 7)) & 0x7e;
  317. #ifdef DEBUG
  318.           if (fpu_debug_level & DL_ARITH)
  319.               printf("__fpu_logn: index to logtbl i=%d(%x)\n", i, i);
  320. #endif
  321.           CPYFPN(&fe->fe_f1, &F);
  322.           /* -F */
  323.           fe->fe_f1.fp_sign = 1;
  324.           /* Y-F */
  325.           d = fpu_add(fe);
  326.           CPYFPN(&fe->fe_f1, d);
  327.  
  328.           /* fe_f2 = 1/F */
  329.           fe->fe_f2.fp_class = FPC_NUM;
  330.           fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[3] = 0;
  331.           fe->fe_f2.fp_exp = logtbl[i].sp_exp;
  332.           fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG));
  333.           fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) |
  334.               (logtbl[i].sp_m1 >> (31 - FP_LG));
  335.           fe->fe_f2.fp_mant[2] = (u_int)(logtbl[i].sp_m1 << (FP_LG + 1));
  336. #ifdef DEBUG
  337.           if (fpu_debug_level & DL_ARITH)
  338.               printf("__fpu_logn: 1/F=(%d,%08x,%08x...)\n", fe->fe_f2.fp_exp,
  339.                        fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
  340. #endif
  341.           /* U = (Y-F) * (1/F) */
  342.           d = fpu_mul(fe);
  343.           CPYFPN(&U, d);
  344.  
  345.           /* KLOG2 = K * ln(2) */
  346.           /* fe_f1 == (fpn)k */
  347.           fpu_explode(fe, &fe->fe_f1, FTYPE_LNG, &k);
  348.           (void)fpu_const(&fe->fe_f2, 0x30 /* ln(2) */);
  349. #ifdef DEBUG
  350.           if (fpu_debug_level & DL_ARITH) {
  351.               printf("__fpu_logn: fp(k)=(%d,%08x,%08x...)\n", fe->fe_f1.fp_exp,
  352.                        fe->fe_f1.fp_mant[0], fe->fe_f1.fp_mant[1]);
  353.               printf("__fpu_logn: ln(2)=(%d,%08x,%08x...)\n", fe->fe_f2.fp_exp,
  354.                        fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
  355.           }
  356. #endif
  357.           /* K * LOGOF2 */
  358.           d = fpu_mul(fe);
  359.           CPYFPN(&KLOG2, d);
  360.  
  361.           /* V=U*U */
  362.           CPYFPN(&fe->fe_f1, &U);
  363.           CPYFPN(&fe->fe_f2, &U);
  364.           d = fpu_mul(fe);
  365.           CPYFPN(&V, d);
  366.  
  367.           /*
  368.            * approximation of LOG(1+U) by
  369.            * (U+V*(A1+V*(A3+V*A5)))+(U*V*(A2+V*(A4+V*A6)))
  370.            */
  371.  
  372.           /* (U+V*(A1+V*(A3+V*A5))) part */
  373.           CPYFPN(&fe->fe_f1, d);
  374.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA5);
  375.           /* V*A5 */
  376.           d = fpu_mul(fe);
  377.  
  378.           CPYFPN(&fe->fe_f1, d);
  379.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA3);
  380.           /* A3+V*A5 */
  381.           d = fpu_add(fe);
  382.  
  383.           CPYFPN(&fe->fe_f1, d);
  384.           CPYFPN(&fe->fe_f2, &V);
  385.           /* V*(A3+V*A5) */
  386.           d = fpu_mul(fe);
  387.  
  388.           CPYFPN(&fe->fe_f1, d);
  389.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA1);
  390.           /* A1+V*(A3+V*A5) */
  391.           d = fpu_add(fe);
  392.  
  393.           CPYFPN(&fe->fe_f1, d);
  394.           CPYFPN(&fe->fe_f2, &V);
  395.           /* V*(A1+V*(A3+V*A5)) */
  396.           d = fpu_mul(fe);
  397.  
  398.           CPYFPN(&fe->fe_f1, d);
  399.           CPYFPN(&fe->fe_f2, &U);
  400.           /* U+V*(A1+V*(A3+V*A5)) */
  401.           d = fpu_add(fe);
  402.  
  403.           CPYFPN(&X, d);
  404.  
  405.           /* (U*V*(A2+V*(A4+V*A6))) part */
  406.           CPYFPN(&fe->fe_f1, &V);
  407.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA6);
  408.           /* V*A6 */
  409.           d = fpu_mul(fe);
  410.           CPYFPN(&fe->fe_f1, d);
  411.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA4);
  412.           /* A4+V*A6 */
  413.           d = fpu_add(fe);
  414.           CPYFPN(&fe->fe_f1, d);
  415.           CPYFPN(&fe->fe_f2, &V);
  416.           /* V*(A4+V*A6) */
  417.           d = fpu_mul(fe);
  418.           CPYFPN(&fe->fe_f1, d);
  419.           fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA2);
  420.           /* A2+V*(A4+V*A6) */
  421.           d = fpu_add(fe);
  422.           CPYFPN(&fe->fe_f1, d);
  423.           CPYFPN(&fe->fe_f2, &V);
  424.           /* V*(A2+V*(A4+V*A6)) */
  425.           d = fpu_mul(fe);
  426.           CPYFPN(&fe->fe_f1, d);
  427.           CPYFPN(&fe->fe_f2, &U);
  428.           /* U*V*(A2+V*(A4+V*A6)) */
  429.           d = fpu_mul(fe);
  430.           CPYFPN(&fe->fe_f1, d);
  431.           i++;
  432.           /* fe_f2 = logtbl[i+1] (== LOG(F)) */
  433.           fe->fe_f2.fp_class = FPC_NUM;
  434.           fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[3] = 0;
  435.           fe->fe_f2.fp_exp = logtbl[i].sp_exp;
  436.           fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG));
  437.           fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) |
  438.               (logtbl[i].sp_m1 >> (31 - FP_LG));
  439.           fe->fe_f2.fp_mant[2] = (logtbl[i].sp_m1 << (FP_LG + 1));
  440. #ifdef DEBUG
  441.           if (fpu_debug_level & DL_ARITH)
  442.               printf("__fpu_logn: ln(F)=(%d,%08x,%08x,...)\n", fe->fe_f2.fp_exp,
  443.                        fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
  444. #endif
  445.           /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
  446.           d = fpu_add(fe);
  447.           CPYFPN(&fe->fe_f1, d);
  448.           CPYFPN(&fe->fe_f2, &X);
  449.           /* LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */
  450.           d = fpu_add(fe);
  451. #ifdef DEBUG
  452.           if (fpu_debug_level & DL_ARITH)
  453.               printf("__fpu_logn: ln(Y)=(%c,%d,%08x,%08x,%08x,%08x)\n",
  454.                        d->fp_sign ? '-' : '+', d->fp_exp,
  455.                        d->fp_mant[0], d->fp_mant[1], d->fp_mant[2], d->fp_mant[3]);
  456. #endif
  457.           CPYFPN(&fe->fe_f1, d);
  458.           CPYFPN(&fe->fe_f2, &KLOG2);
  459.           /* K*LOGOF2+LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */
  460.           d = fpu_add(fe);
  461.     }
  462.  
  463.     return d;
  464. }
  465.  
  466. struct fpn *
  467. fpu_log10(fe)
  468.      struct fpemu *fe;
  469. {
  470.     struct fpn *fp = &fe->fe_f2;
  471.     u_int fpsr;
  472.  
  473.     fpsr = fe->fe_fpsr & ~FPSR_EXCP;    /* clear all exceptions */
  474.  
  475.     if (fp->fp_class >= FPC_NUM) {
  476.           if (fp->fp_sign) {  /* negative number or Inf */
  477.               fp = fpu_newnan(fe);
  478.               fpsr |= FPSR_OPERR;
  479.           } else if (fp->fp_class == FPC_NUM) {
  480.               /* the real work here */
  481.               fp = __fpu_logn(fe);
  482.               if (fp != &fe->fe_f1)
  483.                     CPYFPN(&fe->fe_f1, fp);
  484.               (void)fpu_const(&fe->fe_f2, 0x31 /* ln(10) */);
  485.               fp = fpu_div(fe);
  486.           } /* else if fp == +Inf, return +Inf */
  487.     } else if (fp->fp_class == FPC_ZERO) {
  488.           /* return -Inf */
  489.           fp->fp_class = FPC_INF;
  490.           fp->fp_sign = 1;
  491.           fpsr |= FPSR_DZ;
  492.     } else if (fp->fp_class == FPC_SNAN) {
  493.           fpsr |= FPSR_SNAN;
  494.           fp = fpu_newnan(fe);
  495.     } else {
  496.           fp = fpu_newnan(fe);
  497.     }
  498.  
  499.     fe->fe_fpsr = fpsr;
  500.  
  501.     return fp;
  502. }
  503.  
  504. struct fpn *
  505. fpu_log2(fe)
  506.      struct fpemu *fe;
  507. {
  508.     struct fpn *fp = &fe->fe_f2;
  509.     u_int fpsr;
  510.  
  511.     fpsr = fe->fe_fpsr & ~FPSR_EXCP;    /* clear all exceptions */
  512.  
  513.     if (fp->fp_class >= FPC_NUM) {
  514.           if (fp->fp_sign) {  /* negative number or Inf */
  515.               fp = fpu_newnan(fe);
  516.               fpsr |= FPSR_OPERR;
  517.           } else if (fp->fp_class == FPC_NUM) {
  518.               /* the real work here */
  519.               if (fp->fp_mant[0] == FP_1 && fp->fp_mant[1] == 0 &&
  520.                     fp->fp_mant[2] == 0 && fp->fp_mant[3] == 0) {
  521.                     /* fp == 2.0 ^ exp <--> log2(fp) == exp */
  522.                     fpu_explode(fe, &fe->fe_f3, FTYPE_LNG, &fp->fp_exp);
  523.                     fp = &fe->fe_f3;
  524.               } else {
  525.                     fp = __fpu_logn(fe);
  526.                     if (fp != &fe->fe_f1)
  527.                         CPYFPN(&fe->fe_f1, fp);
  528.                     (void)fpu_const(&fe->fe_f2, 0x30 /* ln(2) */);
  529.                     fp = fpu_div(fe);
  530.               }
  531.           } /* else if fp == +Inf, return +Inf */
  532.     } else if (fp->fp_class == FPC_ZERO) {
  533.           /* return -Inf */
  534.           fp->fp_class = FPC_INF;
  535.           fp->fp_sign = 1;
  536.           fpsr |= FPSR_DZ;
  537.     } else if (fp->fp_class == FPC_SNAN) {
  538.           fpsr |= FPSR_SNAN;
  539.           fp = fpu_newnan(fe);
  540.     } else {
  541.           fp = fpu_newnan(fe);
  542.     }
  543.  
  544.     fe->fe_fpsr = fpsr;
  545.     return fp;
  546. }
  547.  
  548. struct fpn *
  549. fpu_logn(fe)
  550.      struct fpemu *fe;
  551. {
  552.     struct fpn *fp = &fe->fe_f2;
  553.     u_int fpsr;
  554.  
  555.     fpsr = fe->fe_fpsr & ~FPSR_EXCP;    /* clear all exceptions */
  556.  
  557.     if (fp->fp_class >= FPC_NUM) {
  558.           if (fp->fp_sign) {  /* negative number or Inf */
  559.               fp = fpu_newnan(fe);
  560.               fpsr |= FPSR_OPERR;
  561.           } else if (fp->fp_class == FPC_NUM) {
  562.               /* the real work here */
  563.               fp = __fpu_logn(fe);
  564.           } /* else if fp == +Inf, return +Inf */
  565.     } else if (fp->fp_class == FPC_ZERO) {
  566.           /* return -Inf */
  567.           fp->fp_class = FPC_INF;
  568.           fp->fp_sign = 1;
  569.           fpsr |= FPSR_DZ;
  570.     } else if (fp->fp_class == FPC_SNAN) {
  571.           fpsr |= FPSR_SNAN;
  572.           fp = fpu_newnan(fe);
  573.     } else {
  574.           fp = fpu_newnan(fe);
  575.     }
  576.  
  577.     fe->fe_fpsr = fpsr;
  578.  
  579.     return fp;
  580. }
  581.  
  582. struct fpn *
  583. fpu_lognp1(fe)
  584.      struct fpemu *fe;
  585. {
  586.     struct fpn *fp;
  587.  
  588.     /* build a 1.0 */
  589.     fp = fpu_const(&fe->fe_f1, 0x32); /* get 1.0 */
  590.     /* fp = 1.0 + f2 */
  591.     fp = fpu_add(fe);
  592.  
  593.     /* copy the result to the src opr */
  594.     if (&fe->fe_f2 != fp)
  595.           CPYFPN(&fe->fe_f2, fp);
  596.  
  597.     return fpu_logn(fe);
  598. }
  599.